#pragma once
#include <Utilities/glm_ext.h>

namespace PhysLeo {

/**
 * compute the cubic shape derivatives according to different local coordinates.
 * \f$ \frac{\partial{N}}{\partial{X}} \f$ N is the shape function, X is the local coordinates.
 * @param[in] local  the local sampling point, its x y z coordinates must belong to [-1,1]. it is in the local coordinate.
 * @param[out] cubic_shape_derivative  the output cubic shape derivative
 */
template<typename T>
__host__ __device__ void cubicShapeDerivatives(glm::tvec3<T> local, glm::tmat8x3<T>& cubic_shape_derivative)
{
    // shape function
    // N1 = -(X-1)(Y-1)(Z-1)/8
    // N2 = (X+1)(Y-1)(Z-1)/8
    // N3 = -(X+1)(Y+1)(Z-1)/8
    // N4 = (X-1)(Y+1)(Z-1)/8
    // N5 = (X-1)(Y-1)(Z+1)/8
    // N6 = -(X+1)(Y-1)(Z+1)/8
    // N7 = (X+1)(Y+1)(Z+1)/8
    // N8 = -(X-1)(Y+1)(Z+1)/8
    // shape_matrix_(0,i)= \partial{Ni}/\partial{X}

    T X = local[0], Y = local[1], Z = local[2];
    cubic_shape_derivative[0][0] = -(Y - 1)*(Z - 1) / 8;
    cubic_shape_derivative[1][0] = (Y - 1)*(Z - 1) / 8;
    cubic_shape_derivative[2][0] = -(Y + 1)*(Z - 1) / 8;
    cubic_shape_derivative[3][0] = (Y + 1)*(Z - 1) / 8;
    cubic_shape_derivative[4][0] = (Y - 1)*(Z + 1) / 8;
    cubic_shape_derivative[5][0] = -(Y - 1)*(Z + 1) / 8;
    cubic_shape_derivative[6][0] = (Y + 1)*(Z + 1) / 8;
    cubic_shape_derivative[7][0] = -(Y + 1)*(Z + 1) / 8;

    cubic_shape_derivative[0][1] = -(X - 1)*(Z - 1) / 8;
    cubic_shape_derivative[1][1] = (X + 1)*(Z - 1) / 8;
    cubic_shape_derivative[2][1] = -(X + 1)*(Z - 1) / 8;
    cubic_shape_derivative[3][1] = (X - 1)*(Z - 1) / 8;
    cubic_shape_derivative[4][1] = (X - 1)*(Z + 1) / 8;
    cubic_shape_derivative[5][1] = -(X + 1)*(Z + 1) / 8;
    cubic_shape_derivative[6][1] = (X + 1)*(Z + 1) / 8;
    cubic_shape_derivative[7][1] = -(X - 1)*(Z + 1) / 8;

    cubic_shape_derivative[0][2] = -(Y - 1)*(X - 1) / 8;
    cubic_shape_derivative[1][2] = (Y - 1)*(X + 1) / 8;
    cubic_shape_derivative[2][2] = -(Y + 1)*(X + 1) / 8;
    cubic_shape_derivative[3][2] = (Y + 1)*(X - 1) / 8;
    cubic_shape_derivative[4][2] = (Y - 1)*(X - 1) / 8;
    cubic_shape_derivative[5][2] = -(Y - 1)*(X + 1) / 8;
    cubic_shape_derivative[6][2] = (Y + 1)*(X + 1) / 8;
    cubic_shape_derivative[7][2] = -(Y + 1)*(X - 1) / 8;
}

}